Full depth tidal filters from January 1 - February 28, 2018
import qgrid
import xarray as xr
import matplotlib.pyplot as plt
import numpy as np
import scipy.signal as sig
import scipy.interpolate as interp
import pandas as pd
%matplotlib notebook
# Axis January and February
with xr.open_dataset('../Nov11/AxisJanFeb2018.nc') as ds:
print(ds)
grid = pd.DataFrame(uorig)
qgrid.show_grid(grid, show_toolbar=True)
Find specific depth to remove depth values (unreliable data, from visual inspection of initial plots), about 100m at the top
# find specific depth to remove depth values (unreliable data, from visual inspection of initial plots)
# about 100m at the top
def find_nearest(array, value):
array = np.asarray(array)
idx = (np.abs(array - value)).argmin()
return idx # returns index of nearest value
array = ds.depth
upval = 450 # meters
upidx = find_nearest(array, upval)
print("Index at cutoff:", find_nearest(array, upval))
print("Value at cutoff:", abs(array[find_nearest(array, upval)]))
depth = np.array(ds.depth[0:upidx+1]) # remove unwanted upper depths
print("Length of new array: ", len(depth))
Low pass Butterworth filter for 60 hours to remove tides
# low pass Butterworth filter for 60 hours to remove tides
fs = 1 # 1 sample per HOUR for entire time series
fc = 0.0167 # 60 hour low pass filter
Wn = fc / (fs / 2) # normalised cut-off frequencies
b, a = sig.butter(10, Wn, 'lowpass') # digital butterworth filter
w, h = sig.freqz(b, a)
plt.figure(figsize=(8,2)) # check filter response
plt.plot((1/(2*np.pi))*w / 3600, abs(h), label = "60-hour low pass response") # response shown in Hz
plt.xlim(1e-7, 1e-5)
plt.axvline(4.63e-6, color ='red', label = "60-hour cutoff: 4.63e-6 Hz")
plt.xlabel("Frequency [Hz]")
plt.ylabel("Gain")
plt.title("60-hour Butterworth filter frequency response")
plt.legend(loc='best')
plt.show()
Loop to acquire original, filtered, and residual data
# loop to filter tides from all depths
t = len(ds.time) # number of time data points
d = len(depth) # number of depth data points after removing upper portion
days = t/24 # number of days
time = np.linspace(0,days,t) # x-range
uorig = np.empty([t,d]) # empty array for original u data
vorig = np.empty([t,d]) # empty array for original v data
worig = np.empty([t,d]) # empty array for original w data
ulp = np.empty([t,d]) # empty array for low-pass filtered u values
vlp = np.empty([t,d]) # empty array for low-pass filtered v values
wlp = np.empty([t,d]) # empty array for low-pass filtered w values
for j in range(d):
utemp = pd.Series(ds.u[:,j])
uint = utemp.interpolate(method="cubic")
uorig[:,j] = uint # set interpolated data to original array
ulp[:,j] = sig.filtfilt(b, a, uint) # set low pass array values
vtemp = pd.Series(ds.v[:,j])
vint = vtemp.interpolate(method="cubic")
vorig[:,j] = vint # set interpolated data to original array
vlp[:,j] = sig.filtfilt(b, a, vint) # set low pass array values
wtemp = pd.Series(ds.w[:,j])
wint = wtemp.interpolate(method="cubic")
worig[:,j] = wint # set interpolated data to original array
wlp[:,j] = sig.filtfilt(b, a, wint) # set low pass array values
uhp = uorig - ulp
vhp = vorig - vlp
whp = worig - wlp
Plots for u, v, w (original, filtered, and residual) data
# plot u data
fig, (ax1,ax2,ax3) = plt.subplots(3, 1, figsize=(9.5,12), sharex = True)
fig.subplots_adjust(hspace = 0.08)
fig.text(0.5, 0.9, 'Axis ADCP 75 kHz - Jan./Feb. 2018 - u velocities', ha='center', fontsize=12)
fig.text(0.06, 0.5, 'Depth [m]', va='center', rotation='vertical')
fig.text(0.123, 0.887, 'original (NaN interpolated)', va='center')
fig.text(0.123, 0.622, 'filtered', va='center')
fig.text(0.123, 0.359, 'residual', va='center')
fig.text(0.94, 0.5, 'Velocity [m/s]', va='center', rotation='vertical')
fig.text(0.5, 0.08, 'Time [days]', ha='center')
im1 = ax1.pcolormesh(time, -depth, uorig.T, rasterized=True, cmap='RdBu_r', vmin=-0.08, vmax=0.08)
cbar1 = fig.colorbar(im1, ax=ax1, fraction=0.05, pad=0.01, aspect=40)
im2 = ax2.pcolormesh(time, -depth, ulp.T, rasterized=True, cmap='RdBu_r', vmin=-0.08, vmax=0.08)
cbar2 = fig.colorbar(im2, ax=ax2, fraction=0.05, pad=0.01, aspect=40)
im3 = ax3.pcolormesh(time, -depth, uhp.T, rasterized=True, cmap='RdBu_r', vmin=-0.08, vmax=0.08)
cbar3 = fig.colorbar(im3, ax=ax3, fraction=0.05, pad=0.01, aspect=40)
XXX
# plot v data
fig, (ax1,ax2,ax3) = plt.subplots(3, 1, figsize=(9.5,12), sharex = True)
fig.subplots_adjust(hspace = 0.08)
fig.text(0.5, 0.9, 'Axis ADCP 75 kHz - Jan./Feb. 2018 - v velocities', ha='center', fontsize=12)
fig.text(0.06, 0.5, 'Depth [m]', va='center', rotation='vertical')
fig.text(0.123, 0.887, 'original (NaN interpolated)', va='center')
fig.text(0.123, 0.622, 'filtered', va='center')
fig.text(0.123, 0.359, 'residual', va='center')
fig.text(0.94, 0.5, 'Velocity [m/s]', va='center', rotation='vertical')
fig.text(0.5, 0.08, 'Time [days]', ha='center')
im1 = ax1.pcolormesh(time, -depth, vorig.T, rasterized=True, cmap='RdBu_r', vmin=-0.1, vmax=0.1)
cbar1 = fig.colorbar(im1, ax=ax1, fraction=0.05, pad=0.01, aspect=40)
im2 = ax2.pcolormesh(time, -depth, vlp.T, rasterized=True, cmap='RdBu_r', vmin=-0.1, vmax=0.1)
cbar2 = fig.colorbar(im2, ax=ax2, fraction=0.05, pad=0.01, aspect=40)
im3 = ax3.pcolormesh(time, -depth, vhp.T, rasterized=True, cmap='RdBu_r', vmin=-0.1, vmax=0.1)
cbar3 = fig.colorbar(im3, ax=ax3, fraction=0.05, pad=0.01, aspect=40)
XXX
# plot w data
fig, (ax1,ax2,ax3) = plt.subplots(3, 1, figsize=(9.5,12), sharex = True)
fig.subplots_adjust(hspace = 0.08)
fig.text(0.5, 0.9, 'Axis ADCP 75 kHz - Jan./Feb. 2018 - w velocities', ha='center', fontsize=12)
fig.text(0.06, 0.5, 'Depth [m]', va='center', rotation='vertical')
fig.text(0.123, 0.887, 'original (NaN interpolated)', va='center')
fig.text(0.123, 0.622, 'filtered', va='center')
fig.text(0.123, 0.359, 'residual', va='center')
fig.text(0.94, 0.5, 'Velocity [m/s]', va='center', rotation='vertical')
fig.text(0.5, 0.08, 'Time [days]', ha='center')
im1 = ax1.pcolormesh(time, -depth, worig.T, rasterized=True, cmap='RdBu_r', vmin=-0.02, vmax=0.02)
cbar1 = fig.colorbar(im1, ax=ax1, fraction=0.05, pad=0.01, aspect=40)
im2 = ax2.pcolormesh(time, -depth, wlp.T, rasterized=True, cmap='RdBu_r', vmin=-0.02, vmax=0.02)
cbar2 = fig.colorbar(im2, ax=ax2, fraction=0.05, pad=0.01, aspect=40)
im3 = ax3.pcolormesh(time, -depth, whp.T, rasterized=True, cmap='RdBu_r', vmin=-0.02, vmax=0.02)
cbar3 = fig.colorbar(im3, ax=ax3, fraction=0.05, pad=0.01, aspect=40)
XXX
# Axis June and July
with xr.open_dataset('../Nov11/AxisJunJul2018.nc') as ds:
print(ds)
# find specific depth to remove depth values (unreliable data, from visual inspection of initial plots)
# about 100m at the top
def find_nearest(array, value):
array = np.asarray(array)
idx = (np.abs(array - value)).argmin()
return idx # returns index of nearest value
array = ds.depth
upval = 450 # meters
upidx = find_nearest(array, upval)
print("Index at cutoff:", find_nearest(array, upval))
print("Value at cutoff:", abs(array[find_nearest(array, upval)]))
depth = np.array(ds.depth[0:upidx+1]) # remove unwanted upper depths
print("Length of new array: ", len(depth))
# loop to filter tides from all depths
t = len(ds.time) # number of time data points
d = len(depth) # number of depth data points after removing upper portion
days = t/24 # number of days
time = np.linspace(0,days,t) # x-range
uorig = np.empty([t,d]) # empty array for original u data
vorig = np.empty([t,d]) # empty array for original v data
worig = np.empty([t,d]) # empty array for original w data
ulp = np.empty([t,d]) # empty array for low-pass filtered u values
vlp = np.empty([t,d]) # empty array for low-pass filtered v values
wlp = np.empty([t,d]) # empty array for low-pass filtered w values
for j in range(d):
utemp = pd.Series(ds.u[:,j])
uint = utemp.interpolate(method="cubic")
uorig[:,j] = uint # set interpolated data to original array
ulp[:,j] = sig.filtfilt(b, a, uint) # set low pass array values
vtemp = pd.Series(ds.v[:,j])
vint = vtemp.interpolate(method="cubic")
vorig[:,j] = vint # set interpolated data to original array
vlp[:,j] = sig.filtfilt(b, a, vint) # set low pass array values
wtemp = pd.Series(ds.w[:,j])
wint = wtemp.interpolate(method="cubic")
worig[:,j] = wint # set interpolated data to original array
wlp[:,j] = sig.filtfilt(b, a, wint) # set low pass array values
uhp = uorig - ulp
vhp = vorig - vlp
whp = worig - wlp
Plots for u, v, w (original, filtered, and residual) data
# plot u data
fig, (ax1,ax2,ax3) = plt.subplots(3, 1, figsize=(9.5,12), sharex = True)
fig.subplots_adjust(hspace = 0.08)
fig.text(0.5, 0.9, 'Axis ADCP 75 kHz - Jun./Jul. 2018 - u velocities', ha='center', fontsize=12)
fig.text(0.06, 0.5, 'Depth [m]', va='center', rotation='vertical')
fig.text(0.123, 0.887, 'original (NaN interpolated)', va='center')
fig.text(0.123, 0.622, 'filtered', va='center')
fig.text(0.123, 0.359, 'residual', va='center')
fig.text(0.94, 0.5, 'Velocity [m/s]', va='center', rotation='vertical')
fig.text(0.5, 0.08, 'Time [days]', ha='center')
im1 = ax1.pcolormesh(time, -depth, uorig.T, rasterized=True, cmap='RdBu_r', vmin=-0.08, vmax=0.08)
cbar1 = fig.colorbar(im1, ax=ax1, fraction=0.05, pad=0.01, aspect=40)
im2 = ax2.pcolormesh(time, -depth, ulp.T, rasterized=True, cmap='RdBu_r', vmin=-0.08, vmax=0.08)
cbar2 = fig.colorbar(im2, ax=ax2, fraction=0.05, pad=0.01, aspect=40)
im3 = ax3.pcolormesh(time, -depth, uhp.T, rasterized=True, cmap='RdBu_r', vmin=-0.08, vmax=0.08)
cbar3 = fig.colorbar(im3, ax=ax3, fraction=0.05, pad=0.01, aspect=40)
XXX
# plot v data
fig, (ax1,ax2,ax3) = plt.subplots(3, 1, figsize=(9.5,12), sharex = True)
fig.subplots_adjust(hspace = 0.08)
fig.text(0.5, 0.9, 'Axis ADCP 75 kHz - Jun./Jul. 2018 - v velocities', ha='center', fontsize=12)
fig.text(0.06, 0.5, 'Depth [m]', va='center', rotation='vertical')
fig.text(0.123, 0.887, 'original (NaN interpolated)', va='center')
fig.text(0.123, 0.622, 'filtered', va='center')
fig.text(0.123, 0.359, 'residual', va='center')
fig.text(0.94, 0.5, 'Velocity [m/s]', va='center', rotation='vertical')
fig.text(0.5, 0.08, 'Time [days]', ha='center')
im1 = ax1.pcolormesh(time, -depth, vorig.T, rasterized=True, cmap='RdBu_r', vmin=-0.12, vmax=0.12)
cbar1 = fig.colorbar(im1, ax=ax1, fraction=0.05, pad=0.01, aspect=40)
im2 = ax2.pcolormesh(time, -depth, vlp.T, rasterized=True, cmap='RdBu_r', vmin=-0.12, vmax=0.12)
cbar2 = fig.colorbar(im2, ax=ax2, fraction=0.05, pad=0.01, aspect=40)
im3 = ax3.pcolormesh(time, -depth, vhp.T, rasterized=True, cmap='RdBu_r', vmin=-0.12, vmax=0.12)
cbar3 = fig.colorbar(im3, ax=ax3, fraction=0.05, pad=0.01, aspect=40)
XXX
# plot w data
fig, (ax1,ax2,ax3) = plt.subplots(3, 1, figsize=(9.5,12), sharex = True)
fig.subplots_adjust(hspace = 0.08)
fig.text(0.5, 0.9, 'Axis ADCP 75 kHz - Jan./Feb. 2018 - w velocities', ha='center', fontsize=12)
fig.text(0.06, 0.5, 'Depth [m]', va='center', rotation='vertical')
fig.text(0.123, 0.887, 'original (NaN interpolated)', va='center')
fig.text(0.123, 0.622, 'filtered', va='center')
fig.text(0.123, 0.359, 'residual', va='center')
fig.text(0.94, 0.5, 'Velocity [m/s]', va='center', rotation='vertical')
fig.text(0.5, 0.08, 'Time [days]', ha='center')
im1 = ax1.pcolormesh(time, -depth, worig.T, rasterized=True, cmap='RdBu_r', vmin=-0.02, vmax=0.02)
cbar1 = fig.colorbar(im1, ax=ax1, fraction=0.05, pad=0.01, aspect=40)
im2 = ax2.pcolormesh(time, -depth, wlp.T, rasterized=True, cmap='RdBu_r', vmin=-0.02, vmax=0.02)
cbar2 = fig.colorbar(im2, ax=ax2, fraction=0.05, pad=0.01, aspect=40)
im3 = ax3.pcolormesh(time, -depth, whp.T, rasterized=True, cmap='RdBu_r', vmin=-0.02, vmax=0.02)
cbar3 = fig.colorbar(im3, ax=ax3, fraction=0.05, pad=0.01, aspect=40)
XXX